# Loading data and packages, defining functions ----

rm(list = ls())

library(ggplot2)
library(reshape2)
library(stargazer)
library(plm)
library(lmtest)
library(msm)
library(vtable)
library(dplyr)
library(prais)
library(formattable)
library(sparkline)
library(htmlwidgets)
library(webshot)
library(htmltools)
library(corrplot)

# Function for standardizing the data that is used in the panel data models

scaled.data <- function(plm.model){
  temp_data <- scale(plm.model$model)
  temp_data <- data.frame(temp_data)
  temp_data$country <- unlist(lapply(strsplit(rownames(temp_data), "-"), function(x) x[1]))
  temp_data$demyear <- unlist(lapply(strsplit(rownames(temp_data), "-"), function(x) x[2]))
  temp_data <- pdata.frame(temp_data, index=c("country", "demyear"))
  return(temp_data)
}

# Function for calculating the long-term effects of predictors in models that include a lagged dependent variable

lte_se <- function(model, focalvar, lagdv){
  lte <- coef(model)[focalvar] /(1-coef(model)[lagdv])
  tempcov <- vcov(model)[c(focalvar, lagdv), c(focalvar, lagdv)]
  selte <- deltamethod(~ x1/(1-x2), c(coef(model)[focalvar], coef(model)[lagdv]), tempcov, ses = TRUE)
  out <- c(lte-1.96*selte, lte, lte+1.96*selte)
  names(out) <- c("lower", "lte", "upper")
  return(out)
}

# Loading the data, converting to panel data format for the plm package

df <- read.csv("data.csv")
df_original <- df

df <- pdata.frame(df, index=c("country", "demyear"))

# Descriptives table (table 2) ----

desc_table <- st(df, 
   vars = c("closure_fl", "alt_fl", "form_fl", "acc_fl", "pi_z", "epd", "aaph"), 
   labels = c('PSI', "alternation", "formula", "access", 'PI', "established party domination", "party age"),
   out = "return",
   summ = c('mean(x)','sd(x)','min(x)','max(x)'),
   summ.names = c("mean", "SD", "min", "max"))

colnames(desc_table)[1] <- "variable"

vars_for_density <- list(df$closure_fl, df$alt_fl, df$form_fl, df$acc_fl, df$pi_z, df$epd, df$aaph)

desc_table$density <- lapply(vars_for_density, FUN = function(x) as.character(as.tags(sparkline(density(x)$y, type = "line"))))

final_table = as.htmlwidget(formattable(desc_table))

final_table$dependencies <- c(final_table$dependencies, htmlwidgets:::widget_dependencies("sparkline", "sparkline"))

saveWidget(final_table, "bertoa_table_2_descriptives.html", selfcontained = TRUE)

webshot(url = "bertoa_table_2_descriptives.html", file = "bertoa_table_2_descriptives.png", zoom=10)

# Correlations (figure 1) ----

corrmat <- cor(df[, c("closure_fl", "alt_fl", "form_fl", "acc_fl", "pi_z", "epd", "aaph")], use="pairwise.complete.obs")

colnames(corrmat) <- c('PSI', "alternation", "formula", "access", 'PI', "EPD", "party age")
rownames(corrmat) <- c('PSI', "alternation", "formula", "access", 'PI', "EPD", "party age")

png(file = "bertoa_figure_1_correlations.png", width = 6000, height = 4000, bg = "white", pointsize = 160)
corrplot(corrmat, method = "shade", type = "upper", diag = FALSE, addCoef.col = "black", tl.col	= "black", cl.pos	= "n", tl.srt	= 45, col = COL2('RdBu', 200)[-c(1:50, 150:200)])
dev.off()

pdf(file = "bertoa_figure_1_correlations.pdf", width = 6, height = 4, bg = "white")
corrplot(corrmat, method = "shade", type = "upper", diag = FALSE, addCoef.col = "black", tl.col	= "black", cl.pos	= "n", tl.srt	= 45, col = COL2('RdBu', 200)[-c(1:50, 150:200)])
dev.off()

# Time figure (figure 2) ----

df$year5 <- round(df$year/5)*5

temp_out1 <- tapply(df$pi_z[df$year>=1900 & df$year < 2022], df$year5[df$year>=1900 & df$year <= 2022], function(x) mean(x, na.rm=T))

df$closure_z <- scale(as.data.frame(df$closure))

temp_out2 <- tapply(df$closure_z[df$year>=1900 & df$year < 2022], df$year5[df$year>=1900 & df$year < 2022], function(x) mean(x, na.rm=T))

plot_data <- data.frame(year = names(temp_out1), pi_z = temp_out1, closure = temp_out2)

plot_data <- melt(plot_data, id = "year")

levels(plot_data$variable) <- c("party institutionalization", "party system institutionalization")

plot_data$year <- as.numeric(plot_data$year)

ggplot(plot_data, aes(year, value, lty = variable)) +
  geom_line(aes(group = variable)) +
  labs(x = NULL,
       y = "standardized values",
       lty = NULL) +
  theme_bw() +
  theme(panel.grid = element_line(color = "grey95")) +
  theme(legend.position = "bottom", plot.margin = margin(1,1,1,1,"cm"), axis.ticks.y = element_blank())
ggsave("bertoa_figure_2_time.pdf", width = 6, height = 4)

# Granger causality test ----

npercountry <- tapply(1:nrow(df), df$country, length)

df_sub <- df[df$country %in% names(npercountry[which(npercountry > 20)]), ]

test1 <- pgrangertest(closure_fl ~ pi_z, data = df_sub, order = 5, test = "Ztilde")
test1

test2 <- pgrangertest(pi_z ~ closure_fl, data = df_sub, order = 5, test = "Ztilde")
test2

# Overall indices models ----

# Closure

m1_clo_l1 <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag1, data=df, model="within", effect="twoways")
m1_clo_l2 <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag2, data=df, model="within", effect="twoways")
m1_clo_l5 <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag5, data=df, model="within", effect="twoways")

pdwtest(m1_clo_l1)
pdwtest(m1_clo_l2)
pdwtest(m1_clo_l5)

m1_clo_l1_std <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag1, data=scaled.data(m1_clo_l1), model="within", effect="twoways")
m1_clo_l2_std <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag2, data=scaled.data(m1_clo_l2), model="within", effect="twoways")
m1_clo_l5_std <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag5, data=scaled.data(m1_clo_l5), model="within", effect="twoways")

# Party institutionalisation

m1_pi_l1 <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag1, data=df, model="within", effect="twoways")
m1_pi_l2 <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag2, data=df, model="within", effect="twoways")
m1_pi_l5 <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag5, data=df, model="within", effect="twoways")

pdwtest(m1_pi_l1)
pdwtest(m1_pi_l2)
pdwtest(m1_pi_l5)

m1_pi_l1_std <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag1, data=scaled.data(m1_pi_l1), model="within", effect="twoways")
m1_pi_l2_std <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag2, data=scaled.data(m1_pi_l2), model="within", effect="twoways")
m1_pi_l5_std <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag5, data=scaled.data(m1_pi_l5), model="within", effect="twoways")

models_both <- list(m1_clo_l1_std, m1_clo_l2_std, m1_clo_l5_std, m1_pi_l1_std, m1_pi_l2_std, m1_pi_l5_std)

# Table 3

stargazer(models_both,
          type="text",
          se=lapply(models_both, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

# Long term effects (figure 3)

out_clo <- data.frame(rbind(lte_se(m1_clo_l1_std, "pi_z_lag1", "closure_fl_lag1"),
                         lte_se(m1_clo_l2_std, "pi_z_lag2", "closure_fl_lag1"),
                         lte_se(m1_clo_l5_std, "pi_z_lag5", "closure_fl_lag1")), 
                   effect = c("closure", "closure", "closure"),
                   cause = c("PI (t-1)", "PI (t-2)", "PI (t-5)"))

out_pi <- data.frame(rbind(lte_se(m1_pi_l1_std, "closure_fl_lag1", "pi_z_lag1"),
                            lte_se(m1_pi_l2_std, "closure_fl_lag2", "pi_z_lag1"),
                            lte_se(m1_pi_l5_std, "closure_fl_lag5", "pi_z_lag1")), 
                      effect = c("PI", "PI", "PI"),
                      cause = c("closure (t-1)", "closure (t-2)", "closure (t-5)"))

effect_plot_data <- rbind(out_clo, out_pi)

effect_plot_data$effect <- factor(effect_plot_data$effect, levels = c("closure", "PI"), labels = c("closure (t)", "PI (t)"))

effect_plot_data$cause <- factor(effect_plot_data$cause, levels = c("PI (t-1)", "PI (t-2)", "PI (t-5)", "closure (t-1)", "closure (t-2)", "closure (t-5)"))

ggplot(effect_plot_data, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw()
ggsave("bertoa_figure_3_lte.pdf", width=6, height=4)

# Components models ----

## Closure ----

### Alternation ----

m2_alt_epd_l1 <- plm(alt_fl ~ alt_fl_lag1 + epd_lag1, data=df, model="within", effect="twoways")
m2_alt_epd_l2 <- plm(alt_fl ~ alt_fl_lag1 + epd_lag2, data=df, model="within", effect="twoways")
m2_alt_epd_l5 <- plm(alt_fl ~ alt_fl_lag1 + epd_lag5, data=df, model="within", effect="twoways")

m2_alt_aaph_l1 <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag1, data=df, model="within", effect="twoways")
m2_alt_aaph_l2 <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag2, data=df, model="within", effect="twoways")
m2_alt_aaph_l5 <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag5, data=df, model="within", effect="twoways")

pdwtest(m2_alt_epd_l1)
pdwtest(m2_alt_epd_l2)
pdwtest(m2_alt_epd_l5)

pdwtest(m2_alt_aaph_l1)
pdwtest(m2_alt_aaph_l2)
pdwtest(m2_alt_aaph_l5)

m2_alt_epd_l1_std <- plm(alt_fl ~ alt_fl_lag1 + epd_lag1, data=scaled.data(m2_alt_epd_l1), model="within", effect="twoways")
m2_alt_epd_l2_std <- plm(alt_fl ~ alt_fl_lag1 + epd_lag2, data=scaled.data(m2_alt_epd_l2), model="within", effect="twoways")
m2_alt_epd_l5_std <- plm(alt_fl ~ alt_fl_lag1 + epd_lag5, data=scaled.data(m2_alt_epd_l5), model="within", effect="twoways")

m2_alt_aaph_l1_std <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag1, data=scaled.data(m2_alt_aaph_l1), model="within", effect="twoways")
m2_alt_aaph_l2_std <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag2, data=scaled.data(m2_alt_aaph_l2), model="within", effect="twoways")
m2_alt_aaph_l5_std <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag5, data=scaled.data(m2_alt_aaph_l5), model="within", effect="twoways")

table_alt_std <- list(m2_alt_epd_l1_std, m2_alt_epd_l2_std, m2_alt_epd_l5_std, m2_alt_aaph_l1_std, m2_alt_aaph_l2_std, m2_alt_aaph_l5_std)

# Table A1

stargazer(table_alt_std,
          type="text",
          se=lapply(table_alt_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))



out_alt_epd <- data.frame(rbind(lte_se(m2_alt_epd_l1_std, "epd_lag1", "alt_fl_lag1"),
                            lte_se(m2_alt_epd_l2_std, "epd_lag2", "alt_fl_lag1"),
                            lte_se(m2_alt_epd_l5_std, "epd_lag5", "alt_fl_lag1")), 
                      effect = c("alternation", "alternation", "alternation"),
                      cause = c("EPD (t-1)", "EPD (t-2)", "EPD (t-5)"))

out_alt_aaph <- data.frame(rbind(lte_se(m2_alt_aaph_l1_std, "aaph_lag1", "alt_fl_lag1"),
                                 lte_se(m2_alt_aaph_l2_std, "aaph_lag2", "alt_fl_lag1"),
                                 lte_se(m2_alt_aaph_l5_std, "aaph_lag5", "alt_fl_lag1")), 
                           effect = c("alternation", "alternation", "alternation"),
                           cause = c("party age (t-1)", "party age (t-2)", "party age (t-5)"))

### Formula ----

m2_form_epd_l1 <- plm(form_fl ~ form_fl_lag1 + epd_lag1, data=df, model="within", effect="twoways")
m2_form_epd_l2 <- plm(form_fl ~ form_fl_lag1 + epd_lag2, data=df, model="within", effect="twoways")
m2_form_epd_l5 <- plm(form_fl ~ form_fl_lag1 + epd_lag5, data=df, model="within", effect="twoways")

m2_form_aaph_l1 <- plm(form_fl ~ form_fl_lag1 + aaph_lag1, data=df, model="within", effect="twoways")
m2_form_aaph_l2 <- plm(form_fl ~ form_fl_lag1 + aaph_lag2, data=df, model="within", effect="twoways")
m2_form_aaph_l5 <- plm(form_fl ~ form_fl_lag1 + aaph_lag5, data=df, model="within", effect="twoways")

pdwtest(m2_form_epd_l1)
pdwtest(m2_form_epd_l2)
pdwtest(m2_form_epd_l5)

pdwtest(m2_form_aaph_l1)
pdwtest(m2_form_aaph_l2)
pdwtest(m2_form_aaph_l5)

m2_form_epd_l1_std <- plm(form_fl ~ form_fl_lag1 + epd_lag1, data=scaled.data(m2_form_epd_l1), model="within", effect="twoways")
m2_form_epd_l2_std <- plm(form_fl ~ form_fl_lag1 + epd_lag2, data=scaled.data(m2_form_epd_l2), model="within", effect="twoways")
m2_form_epd_l5_std <- plm(form_fl ~ form_fl_lag1 + epd_lag5, data=scaled.data(m2_form_epd_l5), model="within", effect="twoways")

m2_form_aaph_l1_std <- plm(form_fl ~ form_fl_lag1 + aaph_lag1, data=scaled.data(m2_form_aaph_l1), model="within", effect="twoways")
m2_form_aaph_l2_std <- plm(form_fl ~ form_fl_lag1 + aaph_lag2, data=scaled.data(m2_form_aaph_l2), model="within", effect="twoways")
m2_form_aaph_l5_std <- plm(form_fl ~ form_fl_lag1 + aaph_lag5, data=scaled.data(m2_form_aaph_l5), model="within", effect="twoways")

table_form_std <- list(m2_form_epd_l1_std, m2_form_epd_l2_std, m2_form_epd_l5_std, m2_form_aaph_l1_std, m2_form_aaph_l2_std, m2_form_aaph_l5_std)

# Table A2

stargazer(table_form_std,
          type="text",
          se=lapply(table_form_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_form_epd <- data.frame(rbind(lte_se(m2_form_epd_l1_std, "epd_lag1", "form_fl_lag1"),
                                lte_se(m2_form_epd_l2_std, "epd_lag2", "form_fl_lag1"),
                                lte_se(m2_form_epd_l5_std, "epd_lag5", "form_fl_lag1")), 
                          effect = c("formula", "formula", "formula"),
                          cause = c("EPD (t-1)", "EPD (t-2)", "EPD (t-5)"))

out_form_aaph <- data.frame(rbind(lte_se(m2_form_aaph_l1_std, "aaph_lag1", "form_fl_lag1"),
                                 lte_se(m2_form_aaph_l2_std, "aaph_lag2", "form_fl_lag1"),
                                 lte_se(m2_form_aaph_l5_std, "aaph_lag5", "form_fl_lag1")), 
                           effect = c("formula", "formula", "formula"),
                           cause = c("party age (t-1)", "party age (t-2)", "party age (t-5)"))

### Access ----

m2_acc_epd_l1 <- plm(acc_fl ~ acc_fl_lag1 + epd_lag1, data=df, model="within", effect="twoways")
m2_acc_epd_l2 <- plm(acc_fl ~ acc_fl_lag1 + epd_lag2, data=df, model="within", effect="twoways")
m2_acc_epd_l5 <- plm(acc_fl ~ acc_fl_lag1 + epd_lag5, data=df, model="within", effect="twoways")

m2_acc_aaph_l1 <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag1, data=df, model="within", effect="twoways")
m2_acc_aaph_l2 <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag2, data=df, model="within", effect="twoways")
m2_acc_aaph_l5 <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag5, data=df, model="within", effect="twoways")

pdwtest(m2_acc_epd_l1)
pdwtest(m2_acc_epd_l2)
pdwtest(m2_acc_epd_l5)

pdwtest(m2_acc_aaph_l1)
pdwtest(m2_acc_aaph_l2)
pdwtest(m2_acc_aaph_l5)

m2_acc_epd_l1_std <- plm(acc_fl ~ acc_fl_lag1 + epd_lag1, data=scaled.data(m2_acc_epd_l1), model="within", effect="twoways")
m2_acc_epd_l2_std <- plm(acc_fl ~ acc_fl_lag1 + epd_lag2, data=scaled.data(m2_acc_epd_l2), model="within", effect="twoways")
m2_acc_epd_l5_std <- plm(acc_fl ~ acc_fl_lag1 + epd_lag5, data=scaled.data(m2_acc_epd_l5), model="within", effect="twoways")

m2_acc_aaph_l1_std <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag1, data=scaled.data(m2_acc_aaph_l1), model="within", effect="twoways")
m2_acc_aaph_l2_std <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag2, data=scaled.data(m2_acc_aaph_l2), model="within", effect="twoways")
m2_acc_aaph_l5_std <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag5, data=scaled.data(m2_acc_aaph_l5), model="within", effect="twoways")

table_acc_std <- list(m2_acc_epd_l1_std, m2_acc_epd_l2_std, m2_acc_epd_l5_std, m2_acc_aaph_l1_std, m2_acc_aaph_l2_std, m2_acc_aaph_l5_std)

# Table A3

stargazer(table_acc_std,
          type="text",
          se=lapply(table_acc_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_acc_epd <- data.frame(rbind(lte_se(m2_acc_epd_l1_std, "epd_lag1", "acc_fl_lag1"),
                                 lte_se(m2_acc_epd_l2_std, "epd_lag2", "acc_fl_lag1"),
                                 lte_se(m2_acc_epd_l5_std, "epd_lag5", "acc_fl_lag1")), 
                           effect = c("access", "access", "access"),
                           cause = c("EPD (t-1)", "EPD (t-2)", "EPD (t-5)"))

out_acc_aaph <- data.frame(rbind(lte_se(m2_acc_aaph_l1_std, "aaph_lag1", "acc_fl_lag1"),
                                  lte_se(m2_acc_aaph_l2_std, "aaph_lag2", "acc_fl_lag1"),
                                  lte_se(m2_acc_aaph_l5_std, "aaph_lag5", "acc_fl_lag1")), 
                            effect = c("access", "access", "access"),
                            cause = c("party age (t-1)", "party age (t-2)", "party age (t-5)"))

### Overall figure (figure 4) ----

lte_closure_components <- rbind(out_alt_epd, out_alt_aaph, out_form_epd, out_form_aaph, out_acc_epd, out_acc_aaph)

lte_closure_components$effect <- factor(lte_closure_components$effect, levels = c("alternation", "formula", "access"), labels = c(c("alternation (t)", "formula (t)", "access (t)")))

lte_closure_components$cause <- factor(lte_closure_components$cause)

ggplot(lte_closure_components, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("bertoa_figure_4_closure_comp_lte.pdf", width=8, height=4)

## PI ----

### EPD ----

m2_epd_alt_l1 <- plm(epd ~ epd_lag1 + alt_fl_lag1, data=df, model="within", effect="twoways")
m2_epd_alt_l2 <- plm(epd ~ epd_lag1 + alt_fl_lag2, data=df, model="within", effect="twoways")
m2_epd_alt_l5 <- plm(epd ~ epd_lag1 + alt_fl_lag5, data=df, model="within", effect="twoways")

m2_epd_form_l1 <- plm(epd ~ epd_lag1 + form_fl_lag1, data=df, model="within", effect="twoways")
m2_epd_form_l2 <- plm(epd ~ epd_lag1 + form_fl_lag2, data=df, model="within", effect="twoways")
m2_epd_form_l5 <- plm(epd ~ epd_lag1 + form_fl_lag5, data=df, model="within", effect="twoways")

m2_epd_acc_l1 <- plm(epd ~ epd_lag1 + acc_fl_lag1, data=df, model="within", effect="twoways")
m2_epd_acc_l2 <- plm(epd ~ epd_lag1 + acc_fl_lag2, data=df, model="within", effect="twoways")
m2_epd_acc_l5 <- plm(epd ~ epd_lag1 + acc_fl_lag5, data=df, model="within", effect="twoways")

pdwtest(m2_epd_alt_l1)
pdwtest(m2_epd_alt_l2)
pdwtest(m2_epd_alt_l5)

pdwtest(m2_epd_form_l1)
pdwtest(m2_epd_form_l2)
pdwtest(m2_epd_form_l5)

pdwtest(m2_epd_acc_l1)
pdwtest(m2_epd_acc_l2)
pdwtest(m2_epd_acc_l5)

m2_epd_alt_l1_std <- plm(epd ~ epd_lag1 + alt_fl_lag1, data=scaled.data(m2_epd_alt_l1), model="within", effect="twoways")
m2_epd_alt_l2_std <- plm(epd ~ epd_lag1 + alt_fl_lag2, data=scaled.data(m2_epd_alt_l2), model="within", effect="twoways")
m2_epd_alt_l5_std <- plm(epd ~ epd_lag1 + alt_fl_lag5, data=scaled.data(m2_epd_alt_l5), model="within", effect="twoways")

m2_epd_form_l1_std <- plm(epd ~ epd_lag1 + form_fl_lag1, data=scaled.data(m2_epd_form_l1), model="within", effect="twoways")
m2_epd_form_l2_std <- plm(epd ~ epd_lag1 + form_fl_lag2, data=scaled.data(m2_epd_form_l2), model="within", effect="twoways")
m2_epd_form_l5_std <- plm(epd ~ epd_lag1 + form_fl_lag5, data=scaled.data(m2_epd_form_l5), model="within", effect="twoways")

m2_epd_acc_l1_std <- plm(epd ~ epd_lag1 + acc_fl_lag1, data=scaled.data(m2_epd_acc_l1), model="within", effect="twoways")
m2_epd_acc_l2_std <- plm(epd ~ epd_lag1 + acc_fl_lag2, data=scaled.data(m2_epd_acc_l2), model="within", effect="twoways")
m2_epd_acc_l5_std <- plm(epd ~ epd_lag1 + acc_fl_lag5, data=scaled.data(m2_epd_acc_l5), model="within", effect="twoways")

table_epd_std <- list(m2_epd_alt_l1_std, m2_epd_alt_l2_std, m2_epd_alt_l5_std, m2_epd_form_l1_std,  m2_epd_form_l2_std, m2_epd_form_l5_std,  m2_epd_acc_l1_std,  m2_epd_acc_l2_std,  m2_epd_acc_l5_std)

# Table A4

stargazer(table_epd_std,
          type="text",
          se=lapply(table_epd_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_epd_alt <- data.frame(rbind(lte_se(m2_epd_alt_l1_std, "alt_fl_lag1", "epd_lag1"),
                                lte_se(m2_epd_alt_l2_std, "alt_fl_lag2", "epd_lag1"),
                                lte_se(m2_epd_alt_l5_std, "alt_fl_lag5", "epd_lag1")), 
                          effect = c("EPD (t)", "EPD (t)", "EPD (t)"),
                          cause = c("alternation (t-1)", "alternation (t-2)", "alternation (t-5)"))

out_epd_form <- data.frame(rbind(lte_se(m2_epd_form_l1_std, "form_fl_lag1", "epd_lag1"),
                                lte_se(m2_epd_form_l2_std, "form_fl_lag2", "epd_lag1"),
                                lte_se(m2_epd_form_l5_std, "form_fl_lag5", "epd_lag1")), 
                          effect = c("EPD (t)", "EPD (t)", "EPD (t)"),
                          cause = c("formula (t-1)", "formula (t-2)", "formula (t-5)"))

out_epd_acc <- data.frame(rbind(lte_se(m2_epd_acc_l1_std, "acc_fl_lag1", "epd_lag1"),
                                lte_se(m2_epd_acc_l2_std, "acc_fl_lag2", "epd_lag1"),
                                lte_se(m2_epd_acc_l5_std, "acc_fl_lag5", "epd_lag1")), 
                          effect = c("EPD (t)", "EPD (t)", "EPD (t)"),
                          cause = c("access (t-1)", "access (t-2)", "access (t-5)"))

### Party age ----

m2_aaph_alt_l1 <- plm(aaph ~ aaph_lag1 + alt_fl_lag1, data=df, model="within", effect="twoways")
m2_aaph_alt_l2 <- plm(aaph ~ aaph_lag1 + alt_fl_lag2, data=df, model="within", effect="twoways")
m2_aaph_alt_l5 <- plm(aaph ~ aaph_lag1 + alt_fl_lag5, data=df, model="within", effect="twoways")

m2_aaph_form_l1 <- plm(aaph ~ aaph_lag1 + form_fl_lag1, data=df, model="within", effect="twoways")
m2_aaph_form_l2 <- plm(aaph ~ aaph_lag1 + form_fl_lag2, data=df, model="within", effect="twoways")
m2_aaph_form_l5 <- plm(aaph ~ aaph_lag1 + form_fl_lag5, data=df, model="within", effect="twoways")

m2_aaph_acc_l1 <- plm(aaph ~ aaph_lag1 + acc_fl_lag1, data=df, model="within", effect="twoways")
m2_aaph_acc_l2 <- plm(aaph ~ aaph_lag1 + acc_fl_lag2, data=df, model="within", effect="twoways")
m2_aaph_acc_l5 <- plm(aaph ~ aaph_lag1 + acc_fl_lag5, data=df, model="within", effect="twoways")

pdwtest(m2_aaph_alt_l1)
pdwtest(m2_aaph_alt_l2)
pdwtest(m2_aaph_alt_l5)

pdwtest(m2_aaph_form_l1)
pdwtest(m2_aaph_form_l2)
pdwtest(m2_aaph_form_l5)

pdwtest(m2_aaph_acc_l1)
pdwtest(m2_aaph_acc_l2)
pdwtest(m2_aaph_acc_l5)

m2_aaph_alt_l1_std <- plm(aaph ~ aaph_lag1 + alt_fl_lag1, data=scaled.data(m2_aaph_alt_l1), model="within", effect="twoways")
m2_aaph_alt_l2_std <- plm(aaph ~ aaph_lag1 + alt_fl_lag2, data=scaled.data(m2_aaph_alt_l2), model="within", effect="twoways")
m2_aaph_alt_l5_std <- plm(aaph ~ aaph_lag1 + alt_fl_lag5, data=scaled.data(m2_aaph_alt_l5), model="within", effect="twoways")

m2_aaph_form_l1_std <- plm(aaph ~ aaph_lag1 + form_fl_lag1, data=scaled.data(m2_aaph_form_l1), model="within", effect="twoways")
m2_aaph_form_l2_std <- plm(aaph ~ aaph_lag1 + form_fl_lag2, data=scaled.data(m2_aaph_form_l2), model="within", effect="twoways")
m2_aaph_form_l5_std <- plm(aaph ~ aaph_lag1 + form_fl_lag5, data=scaled.data(m2_aaph_form_l5), model="within", effect="twoways")

m2_aaph_acc_l1_std <- plm(aaph ~ aaph_lag1 + acc_fl_lag1, data=scaled.data(m2_aaph_acc_l1), model="within", effect="twoways")
m2_aaph_acc_l2_std <- plm(aaph ~ aaph_lag1 + acc_fl_lag2, data=scaled.data(m2_aaph_acc_l2), model="within", effect="twoways")
m2_aaph_acc_l5_std <- plm(aaph ~ aaph_lag1 + acc_fl_lag5, data=scaled.data(m2_aaph_acc_l5), model="within", effect="twoways")

table_aaph_std <- list(m2_aaph_alt_l1_std, m2_aaph_alt_l2_std, m2_aaph_alt_l5_std, m2_aaph_form_l1_std,  m2_aaph_form_l2_std, m2_aaph_form_l5_std,  m2_aaph_acc_l1_std,  m2_aaph_acc_l2_std,  m2_aaph_acc_l5_std)

# Table A5

stargazer(table_aaph_std,
          type="text",
          se=lapply(table_aaph_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_aaph_alt <- data.frame(rbind(lte_se(m2_aaph_alt_l1_std, "alt_fl_lag1", "aaph_lag1"),
                                lte_se(m2_aaph_alt_l2_std, "alt_fl_lag2", "aaph_lag1"),
                                lte_se(m2_aaph_alt_l5_std, "alt_fl_lag5", "aaph_lag1")), 
                          effect = c("party age (t)", "party age (t)", "party age (t)"),
                          cause = c("alternation (t-1)", "alternation (t-2)", "alternation (t-5)"))

out_aaph_form <- data.frame(rbind(lte_se(m2_aaph_form_l1_std, "form_fl_lag1", "aaph_lag1"),
                                 lte_se(m2_aaph_form_l2_std, "form_fl_lag2", "aaph_lag1"),
                                 lte_se(m2_aaph_form_l5_std, "form_fl_lag5", "aaph_lag1")), 
                           effect = c("party age (t)", "party age (t)", "party age (t)"),
                           cause = c("formula (t-1)", "formula (t-2)", "formula (t-5)"))

out_aaph_acc <- data.frame(rbind(lte_se(m2_aaph_acc_l1_std, "acc_fl_lag1", "aaph_lag1"),
                                lte_se(m2_aaph_acc_l2_std, "acc_fl_lag2", "aaph_lag1"),
                                lte_se(m2_aaph_acc_l5_std, "acc_fl_lag5", "aaph_lag1")), 
                          effect = c("party age (t)", "party age (t)", "party age (t)"),
                          cause = c("access (t-1)", "access (t-2)", "access (t-5)"))


### Overall plot (figure 5) ----

lte_pi_components <- rbind(out_epd_alt, out_epd_form, out_epd_acc, out_aaph_alt, out_aaph_form, out_aaph_acc)

lte_pi_components$effect <- factor(lte_pi_components$effect)

lte_pi_components$cause <- factor(lte_pi_components$cause, levels = c("alternation (t-1)", "alternation (t-2)", "alternation (t-5)", "formula (t-1)", "formula (t-2)", "formula (t-5)", "access (t-1)", "access (t-2)", "access (t-5)")) 

ggplot(lte_pi_components, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("bertoa_figure_5_pi_comp_lte.pdf", width=8, height=4)

# Models with covariates ----

## Overall indices models ----

### Closure ----

m1_clo_l1 <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m1_clo_l2 <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m1_clo_l5 <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

pdwtest(m1_clo_l1)
pdwtest(m1_clo_l2)
pdwtest(m1_clo_l5)

m1_clo_l1_std <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m1_clo_l1), model="within", effect="twoways")
m1_clo_l2_std <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m1_clo_l2), model="within", effect="twoways")
m1_clo_l5_std <- plm(closure_fl ~ closure_fl_lag1 + pi_z_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m1_clo_l5), model="within", effect="twoways")

table_clo_std <- list(m1_clo_l1_std, m1_clo_l2_std, m1_clo_l5_std)

# Table A6

stargazer(table_clo_std,
          type="text",
          se=lapply(table_clo_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

### Party institutionalisation ----

m1_pi_l1 <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m1_pi_l2 <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m1_pi_l5 <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

pdwtest(m1_pi_l1)
pdwtest(m1_pi_l2)
pdwtest(m1_pi_l5)

m1_pi_l1_std <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m1_pi_l1), model="within", effect="twoways")
m1_pi_l2_std <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m1_pi_l2), model="within", effect="twoways")
m1_pi_l5_std <- plm(pi_z ~ pi_z_lag1 + closure_fl_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m1_pi_l5), model="within", effect="twoways")

table_pi_std <- list(m1_pi_l1_std, m1_pi_l2_std, m1_pi_l5_std)

# Table A7

stargazer(table_pi_std,
          type="text",
          se=lapply(table_pi_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

# Long term effects

out_clo <- data.frame(rbind(lte_se(m1_clo_l1_std, "pi_z_lag1", "closure_fl_lag1"),
                            lte_se(m1_clo_l2_std, "pi_z_lag2", "closure_fl_lag1"),
                            lte_se(m1_clo_l5_std, "pi_z_lag5", "closure_fl_lag1")), 
                      effect = c("closure", "closure", "closure"),
                      cause = c("PI (t-1)", "PI (t-2)", "PI (t-5)"))

out_pi <- data.frame(rbind(lte_se(m1_pi_l1_std, "closure_fl_lag1", "pi_z_lag1"),
                           lte_se(m1_pi_l2_std, "closure_fl_lag2", "pi_z_lag1"),
                           lte_se(m1_pi_l5_std, "closure_fl_lag5", "pi_z_lag1")), 
                     effect = c("PI", "PI", "PI"),
                     cause = c("closure (t-1)", "closure (t-2)", "closure (t-5)"))

effect_plot_data <- rbind(out_clo, out_pi)

effect_plot_data$effect <- factor(effect_plot_data$effect, levels = c("closure", "PI"), labels = c("closure (t)", "PI (t)"))

effect_plot_data$cause <- factor(effect_plot_data$cause, levels = c("PI (t-1)", "PI (t-2)", "PI (t-5)", "closure (t-1)", "closure (t-2)", "closure (t-5)"))

# Figure A1

ggplot(effect_plot_data, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw()
ggsave("bertoa_figure_a1_main_lte_covar.pdf")

## Components models ----

### Closure ----

#### Alternation ----

m2_alt_epd_l1 <- plm(alt_fl ~ alt_fl_lag1 + epd_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_alt_epd_l2 <- plm(alt_fl ~ alt_fl_lag1 + epd_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_alt_epd_l5 <- plm(alt_fl ~ alt_fl_lag1 + epd_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

m2_alt_aaph_l1 <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_alt_aaph_l2 <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_alt_aaph_l5 <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

pdwtest(m2_alt_epd_l1)
pdwtest(m2_alt_epd_l2)
pdwtest(m2_alt_epd_l5)

pdwtest(m2_alt_aaph_l1)
pdwtest(m2_alt_aaph_l2)
pdwtest(m2_alt_aaph_l5)

m2_alt_epd_l1_std <- plm(alt_fl ~ alt_fl_lag1 + epd_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_alt_epd_l1), model="within", effect="twoways")
m2_alt_epd_l2_std <- plm(alt_fl ~ alt_fl_lag1 + epd_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_alt_epd_l2), model="within", effect="twoways")
m2_alt_epd_l5_std <- plm(alt_fl ~ alt_fl_lag1 + epd_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_alt_epd_l5), model="within", effect="twoways")

m2_alt_aaph_l1_std <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_alt_aaph_l1), model="within", effect="twoways")
m2_alt_aaph_l2_std <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_alt_aaph_l2), model="within", effect="twoways")
m2_alt_aaph_l5_std <- plm(alt_fl ~ alt_fl_lag1 + aaph_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_alt_aaph_l5), model="within", effect="twoways")

table_alt_std <- list(m2_alt_epd_l1_std, m2_alt_epd_l2_std, m2_alt_epd_l5_std, m2_alt_aaph_l1_std, m2_alt_aaph_l2_std, m2_alt_aaph_l5_std)

# Table A8

stargazer(table_alt_std,
          type="text",
          se=lapply(table_alt_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_alt_epd <- data.frame(rbind(lte_se(m2_alt_epd_l1_std, "epd_lag1", "alt_fl_lag1"),
                                lte_se(m2_alt_epd_l2_std, "epd_lag2", "alt_fl_lag1"),
                                lte_se(m2_alt_epd_l5_std, "epd_lag5", "alt_fl_lag1")), 
                          effect = c("alternation", "alternation", "alternation"),
                          cause = c("EPD (t-1)", "EPD (t-2)", "EPD (t-5)"))

out_alt_aaph <- data.frame(rbind(lte_se(m2_alt_aaph_l1_std, "aaph_lag1", "alt_fl_lag1"),
                                 lte_se(m2_alt_aaph_l2_std, "aaph_lag2", "alt_fl_lag1"),
                                 lte_se(m2_alt_aaph_l5_std, "aaph_lag5", "alt_fl_lag1")), 
                           effect = c("alternation", "alternation", "alternation"),
                           cause = c("party age (t-1)", "party age (t-2)", "party age (t-5)"))

#### Formula ----

m2_form_epd_l1 <- plm(form_fl ~ form_fl_lag1 + epd_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_form_epd_l2 <- plm(form_fl ~ form_fl_lag1 + epd_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_form_epd_l5 <- plm(form_fl ~ form_fl_lag1 + epd_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

m2_form_aaph_l1 <- plm(form_fl ~ form_fl_lag1 + aaph_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_form_aaph_l2 <- plm(form_fl ~ form_fl_lag1 + aaph_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_form_aaph_l5 <- plm(form_fl ~ form_fl_lag1 + aaph_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

pdwtest(m2_form_epd_l1)
pdwtest(m2_form_epd_l2)
pdwtest(m2_form_epd_l5)

pdwtest(m2_form_aaph_l1)
pdwtest(m2_form_aaph_l2)
pdwtest(m2_form_aaph_l5)

m2_form_epd_l1_std <- plm(form_fl ~ form_fl_lag1 + epd_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_form_epd_l1), model="within", effect="twoways")
m2_form_epd_l2_std <- plm(form_fl ~ form_fl_lag1 + epd_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_form_epd_l2), model="within", effect="twoways")
m2_form_epd_l5_std <- plm(form_fl ~ form_fl_lag1 + epd_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_form_epd_l5), model="within", effect="twoways")

m2_form_aaph_l1_std <- plm(form_fl ~ form_fl_lag1 + aaph_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_form_aaph_l1), model="within", effect="twoways")
m2_form_aaph_l2_std <- plm(form_fl ~ form_fl_lag1 + aaph_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_form_aaph_l2), model="within", effect="twoways")
m2_form_aaph_l5_std <- plm(form_fl ~ form_fl_lag1 + aaph_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_form_aaph_l5), model="within", effect="twoways")

table_form_std <- list(m2_form_epd_l1_std, m2_form_epd_l2_std, m2_form_epd_l5_std, m2_form_aaph_l1_std, m2_form_aaph_l2_std, m2_form_aaph_l5_std)

# Table A9

stargazer(table_form_std,
          type="text",
          se=lapply(table_form_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_form_epd <- data.frame(rbind(lte_se(m2_form_epd_l1_std, "epd_lag1", "form_fl_lag1"),
                                 lte_se(m2_form_epd_l2_std, "epd_lag2", "form_fl_lag1"),
                                 lte_se(m2_form_epd_l5_std, "epd_lag5", "form_fl_lag1")), 
                           effect = c("formula", "formula", "formula"),
                           cause = c("EPD (t-1)", "EPD (t-2)", "EPD (t-5)"))

out_form_aaph <- data.frame(rbind(lte_se(m2_form_aaph_l1_std, "aaph_lag1", "form_fl_lag1"),
                                  lte_se(m2_form_aaph_l2_std, "aaph_lag2", "form_fl_lag1"),
                                  lte_se(m2_form_aaph_l5_std, "aaph_lag5", "form_fl_lag1")), 
                            effect = c("formula", "formula", "formula"),
                            cause = c("party age (t-1)", "party age (t-2)", "party age (t-5)"))

#### Access ----

m2_acc_epd_l1 <- plm(acc_fl ~ acc_fl_lag1 + epd_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_acc_epd_l2 <- plm(acc_fl ~ acc_fl_lag1 + epd_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_acc_epd_l5 <- plm(acc_fl ~ acc_fl_lag1 + epd_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

m2_acc_aaph_l1 <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_acc_aaph_l2 <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_acc_aaph_l5 <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

pdwtest(m2_acc_epd_l1)
pdwtest(m2_acc_epd_l2)
pdwtest(m2_acc_epd_l5)

pdwtest(m2_acc_aaph_l1)
pdwtest(m2_acc_aaph_l2)
pdwtest(m2_acc_aaph_l5)

m2_acc_epd_l1_std <- plm(acc_fl ~ acc_fl_lag1 + epd_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_acc_epd_l1), model="within", effect="twoways")
m2_acc_epd_l2_std <- plm(acc_fl ~ acc_fl_lag1 + epd_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_acc_epd_l2), model="within", effect="twoways")
m2_acc_epd_l5_std <- plm(acc_fl ~ acc_fl_lag1 + epd_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_acc_epd_l5), model="within", effect="twoways")

m2_acc_aaph_l1_std <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_acc_aaph_l1), model="within", effect="twoways")
m2_acc_aaph_l2_std <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_acc_aaph_l2), model="within", effect="twoways")
m2_acc_aaph_l5_std <- plm(acc_fl ~ acc_fl_lag1 + aaph_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_acc_aaph_l5), model="within", effect="twoways")

table_acc_std <- list(m2_acc_epd_l1_std, m2_acc_epd_l2_std, m2_acc_epd_l5_std, m2_acc_aaph_l1_std, m2_acc_aaph_l2_std, m2_acc_aaph_l5_std)

# Table A10

stargazer(table_acc_std,
          type="text",
          se=lapply(table_acc_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_acc_epd <- data.frame(rbind(lte_se(m2_acc_epd_l1_std, "epd_lag1", "acc_fl_lag1"),
                                lte_se(m2_acc_epd_l2_std, "epd_lag2", "acc_fl_lag1"),
                                lte_se(m2_acc_epd_l5_std, "epd_lag5", "acc_fl_lag1")), 
                          effect = c("access", "access", "access"),
                          cause = c("EPD (t-1)", "EPD (t-2)", "EPD (t-5)"))

out_acc_aaph <- data.frame(rbind(lte_se(m2_acc_aaph_l1_std, "aaph_lag1", "acc_fl_lag1"),
                                 lte_se(m2_acc_aaph_l2_std, "aaph_lag2", "acc_fl_lag1"),
                                 lte_se(m2_acc_aaph_l5_std, "aaph_lag5", "acc_fl_lag1")), 
                           effect = c("access", "access", "access"),
                           cause = c("party age (t-1)", "party age (t-2)", "party age (t-5)"))

#### Overall plot (figure A2) ----

lte_closure_components <- rbind(out_alt_epd, out_alt_aaph, out_form_epd, out_form_aaph, out_acc_epd, out_acc_aaph)

lte_closure_components$effect <- factor(lte_closure_components$effect, levels = c("alternation", "formula", "access"), labels = c(c("alternation (t)", "formula (t)", "access (t)")))

lte_closure_components$cause <- factor(lte_closure_components$cause)

ggplot(lte_closure_components, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("bertoa_figure_a2_closure_comp_lte_covar.pdf", width=8, height=4)


### PI ----

#### EPD ----

m2_epd_alt_l1 <- plm(epd ~ epd_lag1 + alt_fl_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_epd_alt_l2 <- plm(epd ~ epd_lag1 + alt_fl_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_epd_alt_l5 <- plm(epd ~ epd_lag1 + alt_fl_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

m2_epd_form_l1 <- plm(epd ~ epd_lag1 + form_fl_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_epd_form_l2 <- plm(epd ~ epd_lag1 + form_fl_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_epd_form_l5 <- plm(epd ~ epd_lag1 + form_fl_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

m2_epd_acc_l1 <- plm(epd ~ epd_lag1 + acc_fl_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_epd_acc_l2 <- plm(epd ~ epd_lag1 + acc_fl_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_epd_acc_l5 <- plm(epd ~ epd_lag1 + acc_fl_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

pdwtest(m2_epd_alt_l1)
pdwtest(m2_epd_alt_l2)
pdwtest(m2_epd_alt_l5)

pdwtest(m2_epd_form_l1)
pdwtest(m2_epd_form_l2)
pdwtest(m2_epd_form_l5)

pdwtest(m2_epd_acc_l1)
pdwtest(m2_epd_acc_l2)
pdwtest(m2_epd_acc_l5)

m2_epd_alt_l1_std <- plm(epd ~ epd_lag1 + alt_fl_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_alt_l1), model="within", effect="twoways")
m2_epd_alt_l2_std <- plm(epd ~ epd_lag1 + alt_fl_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_alt_l2), model="within", effect="twoways")
m2_epd_alt_l5_std <- plm(epd ~ epd_lag1 + alt_fl_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_alt_l5), model="within", effect="twoways")

m2_epd_form_l1_std <- plm(epd ~ epd_lag1 + form_fl_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_form_l1), model="within", effect="twoways")
m2_epd_form_l2_std <- plm(epd ~ epd_lag1 + form_fl_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_form_l2), model="within", effect="twoways")
m2_epd_form_l5_std <- plm(epd ~ epd_lag1 + form_fl_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_form_l5), model="within", effect="twoways")

m2_epd_acc_l1_std <- plm(epd ~ epd_lag1 + acc_fl_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_acc_l1), model="within", effect="twoways")
m2_epd_acc_l2_std <- plm(epd ~ epd_lag1 + acc_fl_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_acc_l2), model="within", effect="twoways")
m2_epd_acc_l5_std <- plm(epd ~ epd_lag1 + acc_fl_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_epd_acc_l5), model="within", effect="twoways")

table_epd_std <- list(m2_epd_alt_l1_std, m2_epd_alt_l2_std, m2_epd_alt_l5_std, m2_epd_form_l1_std,  m2_epd_form_l2_std, m2_epd_form_l5_std,  m2_epd_acc_l1_std,  m2_epd_acc_l2_std,  m2_epd_acc_l5_std)

# Table A11

stargazer(table_epd_std,
          type="text",
          se=lapply(table_epd_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_epd_alt <- data.frame(rbind(lte_se(m2_epd_alt_l1_std, "alt_fl_lag1", "epd_lag1"),
                                lte_se(m2_epd_alt_l2_std, "alt_fl_lag2", "epd_lag1"),
                                lte_se(m2_epd_alt_l5_std, "alt_fl_lag5", "epd_lag1")), 
                          effect = c("EPD (t)", "EPD (t)", "EPD (t)"),
                          cause = c("alternation (t-1)", "alternation (t-2)", "alternation (t-5)"))

out_epd_form <- data.frame(rbind(lte_se(m2_epd_form_l1_std, "form_fl_lag1", "epd_lag1"),
                                 lte_se(m2_epd_form_l2_std, "form_fl_lag2", "epd_lag1"),
                                 lte_se(m2_epd_form_l5_std, "form_fl_lag5", "epd_lag1")), 
                           effect = c("EPD (t)", "EPD (t)", "EPD (t)"),
                           cause = c("formula (t-1)", "formula (t-2)", "formula (t-5)"))

out_epd_acc <- data.frame(rbind(lte_se(m2_epd_acc_l1_std, "acc_fl_lag1", "epd_lag1"),
                                lte_se(m2_epd_acc_l2_std, "acc_fl_lag2", "epd_lag1"),
                                lte_se(m2_epd_acc_l5_std, "acc_fl_lag5", "epd_lag1")), 
                          effect = c("EPD (t)", "EPD (t)", "EPD (t)"),
                          cause = c("access (t-1)", "access (t-2)", "access (t-5)"))

#### Party age ----

m2_aaph_alt_l1 <- plm(aaph ~ aaph_lag1 + alt_fl_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_aaph_alt_l2 <- plm(aaph ~ aaph_lag1 + alt_fl_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_aaph_alt_l5 <- plm(aaph ~ aaph_lag1 + alt_fl_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

m2_aaph_form_l1 <- plm(aaph ~ aaph_lag1 + form_fl_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_aaph_form_l2 <- plm(aaph ~ aaph_lag1 + form_fl_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_aaph_form_l5 <- plm(aaph ~ aaph_lag1 + form_fl_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

m2_aaph_acc_l1 <- plm(aaph ~ aaph_lag1 + acc_fl_lag1 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_aaph_acc_l2 <- plm(aaph ~ aaph_lag1 + acc_fl_lag2 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")
m2_aaph_acc_l5 <- plm(aaph ~ aaph_lag1 + acc_fl_lag5 + lsq + log10_gdp + enpp + pol, data=df, model="within", effect="twoways")

pdwtest(m2_aaph_alt_l1)
pdwtest(m2_aaph_alt_l2)
pdwtest(m2_aaph_alt_l5)

pdwtest(m2_aaph_form_l1)
pdwtest(m2_aaph_form_l2)
pdwtest(m2_aaph_form_l5)

pdwtest(m2_aaph_acc_l1)
pdwtest(m2_aaph_acc_l2)
pdwtest(m2_aaph_acc_l5)

m2_aaph_alt_l1_std <- plm(aaph ~ aaph_lag1 + alt_fl_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_alt_l1), model="within", effect="twoways")
m2_aaph_alt_l2_std <- plm(aaph ~ aaph_lag1 + alt_fl_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_alt_l2), model="within", effect="twoways")
m2_aaph_alt_l5_std <- plm(aaph ~ aaph_lag1 + alt_fl_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_alt_l5), model="within", effect="twoways")

m2_aaph_form_l1_std <- plm(aaph ~ aaph_lag1 + form_fl_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_form_l1), model="within", effect="twoways")
m2_aaph_form_l2_std <- plm(aaph ~ aaph_lag1 + form_fl_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_form_l2), model="within", effect="twoways")
m2_aaph_form_l5_std <- plm(aaph ~ aaph_lag1 + form_fl_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_form_l5), model="within", effect="twoways")

m2_aaph_acc_l1_std <- plm(aaph ~ aaph_lag1 + acc_fl_lag1 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_acc_l1), model="within", effect="twoways")
m2_aaph_acc_l2_std <- plm(aaph ~ aaph_lag1 + acc_fl_lag2 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_acc_l2), model="within", effect="twoways")
m2_aaph_acc_l5_std <- plm(aaph ~ aaph_lag1 + acc_fl_lag5 + lsq + log10_gdp + enpp + pol, data=scaled.data(m2_aaph_acc_l5), model="within", effect="twoways")

table_aaph_std <- list(m2_aaph_alt_l1_std, m2_aaph_alt_l2_std, m2_aaph_alt_l5_std, m2_aaph_form_l1_std,  m2_aaph_form_l2_std, m2_aaph_form_l5_std,  m2_aaph_acc_l1_std,  m2_aaph_acc_l2_std,  m2_aaph_acc_l5_std)

# Table A12

stargazer(table_aaph_std,
          type="text",
          se=lapply(table_aaph_std, function(x) coeftest(x, plm::vcovBK(x, type="HC0"))[,2]),
          keep.stat = c("n", "adj.rsq"),
          star.cutoffs = c(0.05, 0.01, 0.001))

out_aaph_alt <- data.frame(rbind(lte_se(m2_aaph_alt_l1_std, "alt_fl_lag1", "aaph_lag1"),
                                 lte_se(m2_aaph_alt_l2_std, "alt_fl_lag2", "aaph_lag1"),
                                 lte_se(m2_aaph_alt_l5_std, "alt_fl_lag5", "aaph_lag1")), 
                           effect = c("party age (t)", "party age (t)", "party age (t)"),
                           cause = c("alternation (t-1)", "alternation (t-2)", "alternation (t-5)"))

out_aaph_form <- data.frame(rbind(lte_se(m2_aaph_form_l1_std, "form_fl_lag1", "aaph_lag1"),
                                  lte_se(m2_aaph_form_l2_std, "form_fl_lag2", "aaph_lag1"),
                                  lte_se(m2_aaph_form_l5_std, "form_fl_lag5", "aaph_lag1")), 
                            effect = c("party age (t)", "party age (t)", "party age (t)"),
                            cause = c("formula (t-1)", "formula (t-2)", "formula (t-5)"))

out_aaph_acc <- data.frame(rbind(lte_se(m2_aaph_acc_l1_std, "acc_fl_lag1", "aaph_lag1"),
                                 lte_se(m2_aaph_acc_l2_std, "acc_fl_lag2", "aaph_lag1"),
                                 lte_se(m2_aaph_acc_l5_std, "acc_fl_lag5", "aaph_lag1")), 
                           effect = c("party age (t)", "party age (t)", "party age (t)"),
                           cause = c("access (t-1)", "access (t-2)", "access (t-5)"))


#### Overall plot (figure A3) ----

lte_pi_components <- rbind(out_epd_alt, out_epd_form, out_epd_acc, out_aaph_alt, out_aaph_form, out_aaph_acc)

lte_pi_components$effect <- factor(lte_pi_components$effect)

lte_pi_components$cause <- factor(lte_pi_components$cause, levels = c("alternation (t-1)", "alternation (t-2)", "alternation (t-5)", "formula (t-1)", "formula (t-2)", "formula (t-5)", "access (t-1)", "access (t-2)", "access (t-5)")) 

ggplot(lte_pi_components, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("bertoa_figure_a3_pi_comp_lte_covar.pdf", width=8, height=4)


# PanelAR ----

df <- df_original

# Standardize the data

df[, c("closure_fl", "pi_z", "closure_fl_lag1", "pi_z_lag1", "closure_fl_lag2", "pi_z_lag2", "closure_fl_lag5", "pi_z_lag5")] <- scale(df[, c("closure_fl", "pi_z", "closure_fl_lag1", "pi_z_lag1", "closure_fl_lag2", "pi_z_lag2", "closure_fl_lag5", "pi_z_lag5")])


m_par_clo_pi_l1 <- prais_winsten(closure_fl ~ pi_z_lag1 + country + as.factor(demyear), data = df, index = c("country", "demyear"))

m_par_clo_pi_l2 <- prais_winsten(closure_fl ~ pi_z_lag2 + country + as.factor(demyear), data = df[complete.cases(df$pi_z_lag2), ], index = c("country", "demyear"))

m_par_clo_pi_l5 <- prais_winsten(closure_fl ~ pi_z_lag5 + country + as.factor(demyear), data = df[complete.cases(df$pi_z_lag5), ], index = c("country", "demyear"))

m_par_pi_clo_l1 <- prais_winsten(pi_z ~ closure_fl_lag1 + country + as.factor(demyear), data = df, index = c("country", "demyear"))

m_par_pi_clo_l2 <- prais_winsten(pi_z ~ closure_fl_lag2 + country + as.factor(demyear), data = df[complete.cases(df$closure_fl_lag2), ], index = c("country", "demyear"))

m_par_pi_clo_l5 <- prais_winsten(pi_z ~ closure_fl_lag5 + country + as.factor(demyear), data = df[complete.cases(df$closure_fl_lag5), ], index = c("country", "demyear"))

panel_ar_justindices <- list(m_par_clo_pi_l1, m_par_clo_pi_l2, m_par_clo_pi_l5, m_par_pi_clo_l1, m_par_pi_clo_l2, m_par_pi_clo_l5)

# Durbin watson statistic

lapply(panel_ar_justindices, function(x) summary(x)$dw[2])

models_out <- lapply(panel_ar_justindices, function(x) coeftest(x, vcov. = vcovPC(x, pairwise = TRUE))[1:2, 1:4])

# Results for table 4

# Intercepts
round(do.call(rbind, lapply(models_out, function(x) x[1,]))[, 1], 3)
round(do.call(rbind, lapply(models_out, function(x) x[1,]))[, 2], 3)

# Coefficients

round(do.call(rbind, lapply(models_out, function(x) x[2,]))[, 1], 3)
round(do.call(rbind, lapply(models_out, function(x) x[2,]))[, 2], 3)

# Number of observations

lapply(panel_ar_justindices, function(x) nrow(x$model))

# Adjusted R-squared

round(unlist(lapply(panel_ar_justindices, function(x) summary(x)$adj.r.squared)), 2)

# Rho

round(unlist(lapply(panel_ar_justindices, function(x) last(summary(x)$rho))), 2)

# Alternative measures for PI and PSI ----

df <- pdata.frame(df_original, index=c("country", "demyear"))

## Overall correlations ----

cor.test(df$closure_fl, df$enpp)
cor.test(df$closure_fl, df$tev)

cor.test(df$pi_z, df$v2xps_party)
cor.test(df$pi_z, df$v2psorgs)
cor.test(df$pi_z, df$v2psprlnks)

cor.test(df$aaph, df$v2psorgs)

## ENPP for closure ----

# Closure

m1_clo_l1 <- plm(enpp ~ enpp_lag1 + pi_z_lag1, data=df, model="within", effect="twoways")
m1_clo_l2 <- plm(enpp ~ enpp_lag1 + pi_z_lag2, data=df, model="within", effect="twoways")
m1_clo_l5 <- plm(enpp ~ enpp_lag1 + pi_z_lag5, data=df, model="within", effect="twoways")

m1_clo_l1_std <- plm(enpp ~ enpp_lag1 + pi_z_lag1, data=scaled.data(m1_clo_l1), model="within", effect="twoways")
m1_clo_l2_std <- plm(enpp ~ enpp_lag1 + pi_z_lag2, data=scaled.data(m1_clo_l2), model="within", effect="twoways")
m1_clo_l5_std <- plm(enpp ~ enpp_lag1 + pi_z_lag5, data=scaled.data(m1_clo_l5), model="within", effect="twoways")

# Party institutionalisation

m1_pi_l1 <- plm(pi_z ~ pi_z_lag1 + enpp_lag1, data=df, model="within", effect="twoways")
m1_pi_l2 <- plm(pi_z ~ pi_z_lag1 + enpp_lag2, data=df, model="within", effect="twoways")
m1_pi_l5 <- plm(pi_z ~ pi_z_lag1 + enpp_lag5, data=df, model="within", effect="twoways")

m1_pi_l1_std <- plm(pi_z ~ pi_z_lag1 + enpp_lag1, data=scaled.data(m1_pi_l1), model="within", effect="twoways")
m1_pi_l2_std <- plm(pi_z ~ pi_z_lag1 + enpp_lag2, data=scaled.data(m1_pi_l2), model="within", effect="twoways")
m1_pi_l5_std <- plm(pi_z ~ pi_z_lag1 + enpp_lag5, data=scaled.data(m1_pi_l5), model="within", effect="twoways")

# Long term effects

out_clo <- data.frame(rbind(lte_se(m1_clo_l1_std, "pi_z_lag1", "enpp_lag1"),
                            lte_se(m1_clo_l2_std, "pi_z_lag2", "enpp_lag1"),
                            lte_se(m1_clo_l5_std, "pi_z_lag5", "enpp_lag1")), 
                      effect = c("enpp", "enpp", "enpp"),
                      cause = c("PI (t-1)", "PI (t-2)", "PI (t-5)"))

out_pi <- data.frame(rbind(lte_se(m1_pi_l1_std, "enpp_lag1", "pi_z_lag1"),
                           lte_se(m1_pi_l2_std, "enpp_lag2", "pi_z_lag1"),
                           lte_se(m1_pi_l5_std, "enpp_lag5", "pi_z_lag1")), 
                     effect = c("PI", "PI", "PI"),
                     cause = c("ENPP (t-1)", "ENPP (t-2)", "ENPP (t-5)"))

effect_plot_data <- rbind(out_clo, out_pi)

effect_plot_data$effect <- factor(effect_plot_data$effect, levels = c("enpp", "PI"), labels = c("ENPP (t)", "PI (t)"))

effect_plot_data$cause <- factor(effect_plot_data$cause, levels = c("PI (t-1)", "PI (t-2)", "PI (t-5)", "ENPP (t-1)", "ENPP (t-2)", "ENPP (t-5)"))

# Figure A4

ggplot(effect_plot_data, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw()
ggsave("bertoa_figure_a4_main_lte_enpp.pdf", width=6, height=4)

## TEV for closure ----

# Closure

m1_clo_l1 <- plm(tev ~ tev_lag1 + pi_z_lag1, data=df, model="within", effect="twoways")
m1_clo_l2 <- plm(tev ~ tev_lag1 + pi_z_lag2, data=df, model="within", effect="twoways")
m1_clo_l5 <- plm(tev ~ tev_lag1 + pi_z_lag5, data=df, model="within", effect="twoways")

m1_clo_l1_std <- plm(tev ~ tev_lag1 + pi_z_lag1, data=scaled.data(m1_clo_l1), model="within", effect="twoways")
m1_clo_l2_std <- plm(tev ~ tev_lag1 + pi_z_lag2, data=scaled.data(m1_clo_l2), model="within", effect="twoways")
m1_clo_l5_std <- plm(tev ~ tev_lag1 + pi_z_lag5, data=scaled.data(m1_clo_l5), model="within", effect="twoways")

# Party institutionalisation

m1_pi_l1 <- plm(pi_z ~ pi_z_lag1 + tev_lag1, data=df, model="within", effect="twoways")
m1_pi_l2 <- plm(pi_z ~ pi_z_lag1 + tev_lag2, data=df, model="within", effect="twoways")
m1_pi_l5 <- plm(pi_z ~ pi_z_lag1 + tev_lag5, data=df, model="within", effect="twoways")

m1_pi_l1_std <- plm(pi_z ~ pi_z_lag1 + tev_lag1, data=scaled.data(m1_pi_l1), model="within", effect="twoways")
m1_pi_l2_std <- plm(pi_z ~ pi_z_lag1 + tev_lag2, data=scaled.data(m1_pi_l2), model="within", effect="twoways")
m1_pi_l5_std <- plm(pi_z ~ pi_z_lag1 + tev_lag5, data=scaled.data(m1_pi_l5), model="within", effect="twoways")

# Long term effects

out_clo <- data.frame(rbind(lte_se(m1_clo_l1_std, "pi_z_lag1", "tev_lag1"),
                            lte_se(m1_clo_l2_std, "pi_z_lag2", "tev_lag1"),
                            lte_se(m1_clo_l5_std, "pi_z_lag5", "tev_lag1")), 
                      effect = c("tev", "tev", "tev"),
                      cause = c("PI (t-1)", "PI (t-2)", "PI (t-5)"))

out_pi <- data.frame(rbind(lte_se(m1_pi_l1_std, "tev_lag1", "pi_z_lag1"),
                           lte_se(m1_pi_l2_std, "tev_lag2", "pi_z_lag1"),
                           lte_se(m1_pi_l5_std, "tev_lag5", "pi_z_lag1")), 
                     effect = c("PI", "PI", "PI"),
                     cause = c("TEV (t-1)", "TEV (t-2)", "TEV (t-5)"))

effect_plot_data <- rbind(out_clo, out_pi)

effect_plot_data$effect <- factor(effect_plot_data$effect, levels = c("tev", "PI"), labels = c("TEV (t)", "PI (t)"))

effect_plot_data$cause <- factor(effect_plot_data$cause, levels = c("PI (t-1)", "PI (t-2)", "PI (t-5)", "TEV (t-1)", "TEV (t-2)", "TEV (t-5)"))

# Figure A5

ggplot(effect_plot_data, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw()
ggsave("bertoa_figure_a5_main_lte_tev.pdf", width=6, height=4)

## VDEM PI index for PI ----

# Closure

m1_clo_l1 <- plm(closure_fl ~ closure_fl_lag1 + v2xps_party_lag1, data=df, model="within", effect="twoways")
m1_clo_l2 <- plm(closure_fl ~ closure_fl_lag1 + v2xps_party_lag2, data=df, model="within", effect="twoways")
m1_clo_l5 <- plm(closure_fl ~ closure_fl_lag1 + v2xps_party_lag5, data=df, model="within", effect="twoways")

m1_clo_l1_std <- plm(closure_fl ~ closure_fl_lag1 + v2xps_party_lag1, data=scaled.data(m1_clo_l1), model="within", effect="twoways")
m1_clo_l2_std <- plm(closure_fl ~ closure_fl_lag1 + v2xps_party_lag2, data=scaled.data(m1_clo_l2), model="within", effect="twoways")
m1_clo_l5_std <- plm(closure_fl ~ closure_fl_lag1 + v2xps_party_lag5, data=scaled.data(m1_clo_l5), model="within", effect="twoways")

# Party institutionalisation

m1_pi_l1 <- plm(v2xps_party ~ v2xps_party_lag1 + closure_fl_lag1, data=df, model="within", effect="twoways")
m1_pi_l2 <- plm(v2xps_party ~ v2xps_party_lag1 + closure_fl_lag2, data=df, model="within", effect="twoways")
m1_pi_l5 <- plm(v2xps_party ~ v2xps_party_lag1 + closure_fl_lag5, data=df, model="within", effect="twoways")

m1_pi_l1_std <- plm(v2xps_party ~ v2xps_party_lag1 + closure_fl_lag1, data=scaled.data(m1_pi_l1), model="within", effect="twoways")
m1_pi_l2_std <- plm(v2xps_party ~ v2xps_party_lag1 + closure_fl_lag2, data=scaled.data(m1_pi_l2), model="within", effect="twoways")
m1_pi_l5_std <- plm(v2xps_party ~ v2xps_party_lag1 + closure_fl_lag5, data=scaled.data(m1_pi_l5), model="within", effect="twoways")

# Long term effects

out_clo <- data.frame(rbind(lte_se(m1_clo_l1_std, "v2xps_party_lag1", "closure_fl_lag1"),
                            lte_se(m1_clo_l2_std, "v2xps_party_lag2", "closure_fl_lag1"),
                            lte_se(m1_clo_l5_std, "v2xps_party_lag5", "closure_fl_lag1")), 
                      effect = c("closure", "closure", "closure"),
                      cause = c("VDEM PI (t-1)", "VDEM PI (t-2)", "VDEM PI (t-5)"))

out_pi <- data.frame(rbind(lte_se(m1_pi_l1_std, "closure_fl_lag1", "v2xps_party_lag1"),
                           lte_se(m1_pi_l2_std, "closure_fl_lag2", "v2xps_party_lag1"),
                           lte_se(m1_pi_l5_std, "closure_fl_lag5", "v2xps_party_lag1")), 
                     effect = c("v2xps_party", "v2xps_party", "v2xps_party"),
                     cause = c("closure (t-1)", "closure (t-2)", "closure (t-5)"))

effect_plot_data <- rbind(out_clo, out_pi)

effect_plot_data$effect <- factor(effect_plot_data$effect, levels = c("closure", "v2xps_party"), labels = c("closure (t)", "VDEM PI (t)"))

effect_plot_data$cause <- factor(effect_plot_data$cause, levels = c("VDEM PI (t-1)", "VDEM PI (t-2)", "VDEM PI (t-5)", "closure (t-1)", "closure (t-2)", "closure (t-5)"))

# Figure A6

ggplot(effect_plot_data, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw()
ggsave("bertoa_figure_a6_main_lte_vdemgeneral.pdf", width=6, height=4)

## VDEM party org for PI ----

# Closure

m1_clo_l1 <- plm(closure_fl ~ closure_fl_lag1 + v2psorgs_lag1, data=df, model="within", effect="twoways")
m1_clo_l2 <- plm(closure_fl ~ closure_fl_lag1 + v2psorgs_lag2, data=df, model="within", effect="twoways")
m1_clo_l5 <- plm(closure_fl ~ closure_fl_lag1 + v2psorgs_lag5, data=df, model="within", effect="twoways")

m1_clo_l1_std <- plm(closure_fl ~ closure_fl_lag1 + v2psorgs_lag1, data=scaled.data(m1_clo_l1), model="within", effect="twoways")
m1_clo_l2_std <- plm(closure_fl ~ closure_fl_lag1 + v2psorgs_lag2, data=scaled.data(m1_clo_l2), model="within", effect="twoways")
m1_clo_l5_std <- plm(closure_fl ~ closure_fl_lag1 + v2psorgs_lag5, data=scaled.data(m1_clo_l5), model="within", effect="twoways")

# Party institutionalisation

m1_pi_l1 <- plm(v2psorgs ~ v2psorgs_lag1 + closure_fl_lag1, data=df, model="within", effect="twoways")
m1_pi_l2 <- plm(v2psorgs ~ v2psorgs_lag1 + closure_fl_lag2, data=df, model="within", effect="twoways")
m1_pi_l5 <- plm(v2psorgs ~ v2psorgs_lag1 + closure_fl_lag5, data=df, model="within", effect="twoways")

m1_pi_l1_std <- plm(v2psorgs ~ v2psorgs_lag1 + closure_fl_lag1, data=scaled.data(m1_pi_l1), model="within", effect="twoways")
m1_pi_l2_std <- plm(v2psorgs ~ v2psorgs_lag1 + closure_fl_lag2, data=scaled.data(m1_pi_l2), model="within", effect="twoways")
m1_pi_l5_std <- plm(v2psorgs ~ v2psorgs_lag1 + closure_fl_lag5, data=scaled.data(m1_pi_l5), model="within", effect="twoways")

# Long term effects

out_clo <- data.frame(rbind(lte_se(m1_clo_l1_std, "v2psorgs_lag1", "closure_fl_lag1"),
                            lte_se(m1_clo_l2_std, "v2psorgs_lag2", "closure_fl_lag1"),
                            lte_se(m1_clo_l5_std, "v2psorgs_lag5", "closure_fl_lag1")), 
                      effect = c("closure", "closure", "closure"),
                      cause = c("VDEM PO (t-1)", "VDEM PO (t-2)", "VDEM PO (t-5)"))

out_pi <- data.frame(rbind(lte_se(m1_pi_l1_std, "closure_fl_lag1", "v2psorgs_lag1"),
                           lte_se(m1_pi_l2_std, "closure_fl_lag2", "v2psorgs_lag1"),
                           lte_se(m1_pi_l5_std, "closure_fl_lag5", "v2psorgs_lag1")), 
                     effect = c("v2psorgs", "v2psorgs", "v2psorgs"),
                     cause = c("closure (t-1)", "closure (t-2)", "closure (t-5)"))

effect_plot_data <- rbind(out_clo, out_pi)

effect_plot_data$effect <- factor(effect_plot_data$effect, levels = c("closure", "v2psorgs"), labels = c("closure (t)", "VDEM PO (t)"))

effect_plot_data$cause <- factor(effect_plot_data$cause, levels = c("VDEM PO (t-1)", "VDEM PO (t-2)", "VDEM PO (t-5)", "closure (t-1)", "closure (t-2)", "closure (t-5)"))

# Figure A7

ggplot(effect_plot_data, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw()
ggsave("bertoa_figure_a7_main_lte_vdemporg.pdf", width=6, height=4)

## VDEM linkages for PI ----

# Closure

m1_clo_l1 <- plm(closure_fl ~ closure_fl_lag1 + v2psprlnks_lag1, data=df, model="within", effect="twoways")
m1_clo_l2 <- plm(closure_fl ~ closure_fl_lag1 + v2psprlnks_lag2, data=df, model="within", effect="twoways")
m1_clo_l5 <- plm(closure_fl ~ closure_fl_lag1 + v2psprlnks_lag5, data=df, model="within", effect="twoways")

m1_clo_l1_std <- plm(closure_fl ~ closure_fl_lag1 + v2psprlnks_lag1, data=scaled.data(m1_clo_l1), model="within", effect="twoways")
m1_clo_l2_std <- plm(closure_fl ~ closure_fl_lag1 + v2psprlnks_lag2, data=scaled.data(m1_clo_l2), model="within", effect="twoways")
m1_clo_l5_std <- plm(closure_fl ~ closure_fl_lag1 + v2psprlnks_lag5, data=scaled.data(m1_clo_l5), model="within", effect="twoways")

# Party institutionalisation

m1_pi_l1 <- plm(v2psprlnks ~ v2psprlnks_lag1 + closure_fl_lag1, data=df, model="within", effect="twoways")
m1_pi_l2 <- plm(v2psprlnks ~ v2psprlnks_lag1 + closure_fl_lag2, data=df, model="within", effect="twoways")
m1_pi_l5 <- plm(v2psprlnks ~ v2psprlnks_lag1 + closure_fl_lag5, data=df, model="within", effect="twoways")

m1_pi_l1_std <- plm(v2psprlnks ~ v2psprlnks_lag1 + closure_fl_lag1, data=scaled.data(m1_pi_l1), model="within", effect="twoways")
m1_pi_l2_std <- plm(v2psprlnks ~ v2psprlnks_lag1 + closure_fl_lag2, data=scaled.data(m1_pi_l2), model="within", effect="twoways")
m1_pi_l5_std <- plm(v2psprlnks ~ v2psprlnks_lag1 + closure_fl_lag5, data=scaled.data(m1_pi_l5), model="within", effect="twoways")

# Long term effects

out_clo <- data.frame(rbind(lte_se(m1_clo_l1_std, "v2psprlnks_lag1", "closure_fl_lag1"),
                            lte_se(m1_clo_l2_std, "v2psprlnks_lag2", "closure_fl_lag1"),
                            lte_se(m1_clo_l5_std, "v2psprlnks_lag5", "closure_fl_lag1")), 
                      effect = c("closure", "closure", "closure"),
                      cause = c("VDEM PL (t-1)", "VDEM PL (t-2)", "VDEM PL (t-5)"))

out_pi <- data.frame(rbind(lte_se(m1_pi_l1_std, "closure_fl_lag1", "v2psprlnks_lag1"),
                           lte_se(m1_pi_l2_std, "closure_fl_lag2", "v2psprlnks_lag1"),
                           lte_se(m1_pi_l5_std, "closure_fl_lag5", "v2psprlnks_lag1")), 
                     effect = c("v2psprlnks", "v2psprlnks", "v2psprlnks"),
                     cause = c("closure (t-1)", "closure (t-2)", "closure (t-5)"))

effect_plot_data <- rbind(out_clo, out_pi)

effect_plot_data$effect <- factor(effect_plot_data$effect, levels = c("closure", "v2psprlnks"), labels = c("closure (t)", "VDEM PL (t)"))

effect_plot_data$cause <- factor(effect_plot_data$cause, levels = c("VDEM PL (t-1)", "VDEM PL (t-2)", "VDEM PL (t-5)", "closure (t-1)", "closure (t-2)", "closure (t-5)"))

# Figure A8

ggplot(effect_plot_data, aes(cause, lte, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  geom_errorbar(width = 0.5) +
  facet_wrap(~effect, ncol = 4, scales = "free_x") +
  labs(y = "long-term effect",
       x = "cause",
       title = NULL) +
  theme_bw()
ggsave("bertoa_figure_a8_main_lte_vdemplink.pdf", width=6, height=4)
